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1. Introduction 

Experimental observations have shown that tokamak plasmas rotate spontaneously 
without momentum input [1]. This intrinsic rotation has been the object of recent 
work [H |2] because of its relevance for ITER [3], where the projected momentum input 
from neutral beams is small, and the rotation is expected to be mostly intrinsic. 

The origin of the intrinsic rotation is still unclear. There has been some 
theoretical work in turbulent transport of momentum using gyrokinetic simulations 
[H El |6l [71 m O [ini [m 112], and two main mechanisms have been proposed as candidates 
to explain intrinsic rotation. On the one hand, the momentum pinch due to the Coriolis 
drift [4j has been argued to transport momentum generated in the edge. On the 
other hand, it has been argued that up-down asymmetry generates intrinsic rotation 
[3 IE] . However, neither of these explanations are able to account for all experimental 
observations. The up-down asymmetry is only large in the edge, generating rotation in 
that region that then needs to be transported inwards by the Coriolis pinch. Thus, 
intrinsic rotation in the core could only be explained by the pinch. The pinch of 
momentum is not sufficient because it does not allow the toroidal rotation to change 
sign in the core as is observed experimentally [T3] . 

In this article we present a new model implementable in 5f flux tube simulations 
[m [ISl [ini [IT] . This model is based on the low flow ordering of [18], and self-consistently 
includes higher order contributions. As a result, new drive terms for the intrinsic 
rotation appear that depend on the gradients of the background profiles of density 
and temperature. 

We recast the results from [T8| in a form similar to the equations in the high 
flow ordering [191 EO]- These are the equations that have been implemented in most 
gyrokinetic codes that are employed to study momentum transport. For this reason, 
the new form of the equations is useful to identify the differences with previous models. 
In addition, we discuss how the new contributions drive intrinsic rotation and we show 
that the intrinsic rotation resulting from these new processes depends on density and 
temperature gradients. 

In the remainder of this article we present the model, developed originally in [18] . 
in a form more suitable for 5f flux tube simulation. In Section [2] we give the complete 
model, and in Section [3] we discuss its implications for intrinsic rotation. Appendix [Xj 
contains the details of the transformation from the equations in [18j to the formulation 
in this article. In Appendix [B] we discuss different forms of deriving the radial flux of 
toroidal angular momentum, showing that the flnal form presented here is convenient 
and has many advantages. 

2. Transport of toroidal angular momentum 

The derivation of the transport of toroidal angular momentum in the low flow regime, 
including both turbulence and neoclassical effects, is described in detail in [18]. To 
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simplify the derivation, the extra expansion parameter Bp/B = s/q <^ 1 was employed, 
assuming that the turbulence length scales and amplitudes do not depend strongly on 
Bp/B. Here B is the total magnetic field and Bp its poloidal component, e = r / R is 
the inverse aspect ratio of the flux surface, q is the safety factor, r is the minor radius 
of the flux surface and R is the major radius. The ratio Bp/B is below or around 0.1 
across the core in most tokamaks [e is small near the magnetic axis and q is large near 
the edge). In this section, we review the results of reference [18] and we recast them in 
a more convenient form. 

We assume that the turbulence is electrostatic and that the magnetic fleld is 
axisymmetric, i.e., B = JVC + VC x VV', where ip is the poloidal magnetic flux, ( 
is the toroidal angle, and we use a poloidal angle 6 as our third spatial coordinate. 
With an axisymmetric magnetic fleld, in steady state and in the absence of momentum 

input, the equation that determines the rotation proflle is {{RC- Pi ■'ViP)^)t = 0, where 
•<->■ 

P .= J dP-v' fiMW is the ion stress tensor, M is the ion mass, <^ is the unit vector in 
the toroidal direction, (. . .)^ = (y')~^ J d9dC. (...)/(B ■ V9) is the flux surface average, 
V' = dV/dip = j d6 d( (B ■ V^)~^ is the derivative of the volume with respect to ip, 
and (. . .)t is the coarse grain or "transport" average over the time and length scales of 
the turbulence, assumed much shorter than the transport time scale 5~'^a/vti and the 
minor radius a. Here 5j = p^/a <^ 1 is the ion gyroradius pi over the minor radius of 
the tokamak a, and Vu is the ion thermal speed. Note that we use the prime in v' to 
indicate that the velocity is measured in the laboratory frame. Later we will flnd the 
equations in a convenient rotating frame where the velocity is v = v' — 

In reference [H] we derived a method to calculate ((-RC' Pj ■V'?/')^)t to order 
{B/Bp)5fpiR\Vip\, with pi the ion pressure. We present the method again in a different 
form to make it easier to compare with previous works in the high flow regime [19], [20| . 
In subsection 12. II we explain how we split the distribution function and the electrostatic 
potential into different pieces, and we present the equations to self-consistently obtain 
them. In subsection 12.21 we evaluate ((-RC' P« ■V'?/^)^)t employing the pieces of the 
distribution function and the potential obtained in subsection 12.11 Before presenting all 
the results, we emphasize that our results and order of magnitude estimates are valid for 
6i <^ Bp/B <^ 1, assuming that the turbulence does not scale strongly with Bp/B, and 
for coUisionality in the range 6f <^ qRvii/vti<l [18], where vu is the ion-ion collision 
frequency. 

2.1. Distribution function and electrostatic potential 

The electrostatic potential is composed to the order of interest by the pieces in Table [I] 
[T8] . The axisymmetric long wavelength pieces (f)o{ip,t), 4>i'^{4',9,t) and 02'^ (■?/', 6*, are 
the zeroth, flrst and second order equilibrium pieces of the potential. The lowest order 
component 0o is a flux surface function. The corrections (^5"^ and 02^^ give the electric fleld 
parallel to the flux surface, established to force quasineutrality at long wavelengths (the 
superscript '^'^ refers to neoclassical because these are long wavelength contributions; 
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Table 1. Pieces of the potential: ip 


= '?!>o ^ 


- <i>r + 't>r 




Potential 


Size 


Leng 


;tli scales 


Time scales 




Te/e 


ka 


' 1 


djdt ~ 5fvu/a 




{B/Bp)d,Tje 


fca ~ 


' 1 


d/dt ^ 5fvti/a 




{B/B,)H^Tje 


ka ^ 


' 1 


djdt ^ 5jvti/a 




~ 6,Te/e 


ti^Pi 


- 1 


d/dt ^ vti/a 




4>'^ ~ {BIB^)SlT,le 


k\\R 


- 1 





however, turbulence can affect the final value of 02'^). We need not calculate 02^^ 

because it will not appear in the final expression for {{RC,- Pi ■V'0)^)t- The piece 
0*'^(r,t) is turbulent and includes both axisymmetric components (zonal flow) and non- 
axisymmetric fluctuations. It is small in 5, but it has strong perpendicular gradients, 
i.e., k^pi ~ 1. Its parallel gradient is small, i.e., k\\R ~ 1. The function cf)^^ is calculated 
to order {B/Bp)6fTje, i.e., 0*^ = 0f + 0f with 0f ~ 6iTje and 0*2^ ~ {B/Bp)6fTje. 
It is convenient to keep both pieces together as 0**^ as we do hereafter. 

To write the distribution function it will be useful to consider the reference frame 
that rotates with toroidal angular velocity = —cd^cf)^. In this new reference frame it 
is easier to compare with previous formulations [121 ED]- To shorten the presentation, we 
perform the change of reference frame directly in the gyrokinetic variables. It is possible 
to do so easily because we are expanding in the parameter B/Bp ^ 1 . We first present 
the gyrokinetic variables that we obtained for the laboratory frame and we argue later 
how they must be modified to give the gyrokinetic variables in the rotating frame. In 
[IS] we used as gyrokinetic variables the gyrocenter position R = r + Ri + R2 + . . ., the 
gyrokinetic kinetic energy E = E0 + E1 + E2 + ■ ■ ., the magnetic moment p = fio + + ■ ■ ■ 
and the gyrokinetic gyrophase ip = ipQ + ipi + . . ., where Eq = {v')^/2 is the particle 
kinetic energy in the laboratory frame, po = {v'^Y/2B is the lowest order magnetic 
moment, = arctan(v' ■ e2/v' • ei) is the lowest order gyrophase, Ri = ^^^"v' x b ~ 6ia 
is the first order correction to the gyrocenter position, Ei = Ze{(j) — {(f))i)/M ~ 6iV^^ is 
the first order correction to the gyrokinetic kinetic energy, and the corrections R2 ~ Sfa, 
E2 ~ '^j^'^i, l^i ~ ^i'^ti/^ ~ defined in [21]. Here fli = ZeB/Mc is the ion 

gyrofrequency, ei(r) and 62 (r) are two orthonormal vectors such that ei x 62 = b, and 
(. . .)i = (27r)"^ § d{p{. . .)\■R_^E,^l,t is the gyroaverage holding R, E, p and t fixed. When 
the ion distribution function is written as a function of these gyrokinetic variables, it 
does not depend on the gyrophase ip up to order {Bp/ B)5f{qRi'ii/vti)fMi [HI 121], where 
fMi is the lowest order distribution function that is a Maxwellian. For the magnetic 
moment and the gyrophase, only the first order corrections /ii and ipi are needed because 
the lowest order distribution function does not depend on p or ip. Moreover, in [IS] 
we expand for 1 Bp/B ^ 6i, and the distribution function need only be known to 
order {B/Bp)SffMi- Consequently, the piece of the distribution function that depends 
on the gyrophase, of order {Bp/B)6f{qRh'ii/vu)fMi, is negligible, and the gyrokinetic 
variables R and E only need to be obtained to order {B/ Bp)6fa and {B/Bp)6fv^^, 
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respectively, implying that the corrections R2 and E2 are not needed for the final 
result. To change to the new reference frame, where the velocity is v = v' — RVtf^C,, 
the distribution function that is independent of if has to be written as a function of the 
new gyrokinetic variables R, e and /x, where e is a new variable that will be defined 
shortly. Note that the gyrocenter position and the magnetic moment are the same in 
both reference frames to the order of interest. In the case of /x, the reason is that /i is 
obtained such that its time derivative vanishes, dfi/dt = 0, making its definition unique. 
For R, the reason is that the toroidal rotation has two components, one parallel to 
the magnetic field, Rflt^C ' b = I^(/B ~ {B / Bp)6iVu, and the other perpendicular, 
— hh ■ ^1 = \Vip\fl(^/B ~ 6iVu, and the parallel velocity is larger by B/Bp ^ 1. 
Since the gyrokinetic variable R is to be obtained to order {B/ Bp)6fa, and in R only 
the perpendicular velocity = + — bb • ^) enters to order 6ia, we can safely 

neglect the corrections due to the change of reference frame because they are of order 
Sfa. In contrast, the kinetic energy E as defined in [18] cannot be used in the rotating 
frame because it includes the parallel velocity U|| = v\\ + IQi^/B. We use a new kinetic 
energy variable e that is related to the old kinetic energy variable hj e = E — IQ^^u'/B, 
where u' = ±a/2(-E' — fiB) is the gyrokinetic parallel velocity in the laboratory frame. 

It is easy to check that u = ±^2[6 — fiB + {I /ByQ'^/2] is equal to u = u' — IVtc^/B 
and it is the gyrokinetic parallel velocity in the rotating frame. With this relation, we 
find that another way to interpret the new energy variable 

e=- + ,B^^ (1) 

is realizing that it is the kinetic energy in the rotating frame plus the potential due to 
the centrifugal force. To write expression ([T]) we have used that I / B = R+0[{Bp/ B)'^R\ 
for Bp/B <^ 1. In Appendix |X] we rewrite the results in using the new gyrokinetic 
kinetic energy e. 

The different pieces of the ion distribution function are given in Table [2] [18] . The 
functions /a/j, H'^i-, H^2 axisymmetric long wavelength contributions. The 

Maxwellian 

M 



fMi{^{'R),e) = ni{^{R)) 



3/2 / 



.2irri('i;)(R)). 

is uniform on a fiux surface. The first and second order corrections Hff^ and H^2 ^-'^s 
neoclassical corrections, and they are not the functions F^^ and F^ in [IB] because we 
are now working in the rotating frame. The function H^^ is an axisymmetric piece of 
the distribution function that originates from collisions acting on the ions transported 
by turbulent fiuctuations into the volume between two adjacent fiux surfaces [18]. The 
function fl^ is the turbulent contribution. It will be determined self-consistently up 
to order {B/Bp)6ffM^, i.e., /^"^ = + with /^f ~ and ~ {B/Bp)SUMi. 

It is convenient to combine both pieces of the turbulent distribution function into one 
function fl^. 

The electron distribution function is very similar to the ion distribution function. 
It will have its own gyrokinetic variables that can be easily deduced from the ion 
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Table 2. Pieces of the ion distribution function: + H'^ + + Hf^ + ff". 



Distribution function Size Length scales Time scales 



/Mi(V'(R),^,^) 




fhli 


ka ~ 


1 


d/dt 


~ djvti/a 


i7-(^(R),^^(R), 


e,fi,t) 


{B /Bp)6ifMi 


ka ~ 


1 


d/dt 


~ 5fvti/a 




e,n,t) 


{B/B,)H^f,u 


ka ~ 


1 


d/dt 


~ 5fvti/a 




e,fx,t) 


{B/Bp){vu/qRuii)SffMi 


ka ~ 


1 


d/dt 


~ ^fvti/a 






~ {B/Bp)5^fM^ 


k±Pi 
k\\Rr 


~ 1 
- 1 


d/dt 


~ vti/a 



Table 3. Pieces of the electron distribution function: /e = fue + -^"1 + /e ■ 



Distribution function Size Length scales Time scales 

/Me(V'(R-),e,i) fMe fcfl 1 d/dt^Slvtila 

H,"f(V'(R),0(R),e,Ai,O {B/Bp)5,fMe ka ^ I djdt^Sfvu/a 

/*''(R, e, /i, t) ^ SJmc k^p, - 1 d/dt - wtj/a 



counterparts. To the order of interest in this calculation, the electron distribution 
function is determined by the pieces in Table [31 The long wavelength, axisymmetric 
pieces Jmb and are the lowest order Maxwellian and the first order neoclassical 
correction. The second order long wavelength neoclassical correction is not needed 
for transport of momentum because of the small electron mass. The piece f^^ is the 
short wavelength, turbulent component that will be self-consistently calculated to order 

iB/Bp)5ffMe. 

We now proceed to describe how to find the different pieces of the distribution 
function and the potential. We use the equations in |T8] but we change to the new 
gyrokinetic kinetic energy e. The details of this transformation are contained in 
Appendix |Al 



2.1.1. First order neoclassical distribution function and potential. The equation for 
i/- is 



«b-VR \ + -^fM^ + 



1 dpi ^ 
Pi d^j 



Me 



l_dJ\ 

Tidij 



luf. 



Mi 



where u = ±^2{e — pB + i?2f2^/2) ~ ±A/2(e — fiB) is the gyrokinetic parallel velocity 

and C^P is the linearized ion-ion collision operator. The correction Hf^ gives the 
parallel component of the velocity [22ll23] UiWf^ = b J d^v Hfl'vw = —{cIh/ZeB)d^Pi + 
{knicI'B/Ze{B^)^)dtpTi, where A; is a flux function that depends on the coUisionality and 
the magnetic geometry. 
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The equation for H^^ is similar to (E]) and it is given by 
. / e^r f \ldp (Ms 5\ 1 9TJ lufMe \ 

-Cif {i/,T} - C^{H:n = -'4^uh . (4) 

where m and fig = eB/mc are the electron mass and gyrofrequency, E"^ is the electric 
field driven by the transformer, Ce^ is the linearized electron-electron collision operator 
and Cgf is the linearized electron-ion collision operator. The lowest order solution for 
ifgf is the Maxwell-Boltzmann response (e05"^/Te)/Me ~ {B / Bp)5ifMe- The rest of 
the terms are small because they are of order {B / Bp)5efMe ~ {B / Bp)^Jm/ M5ifMi ^ 
{B/ Bp)SifMe, where 6e = Pe/a is the ratio between the electron gyroradius pe and the 
minor radius a. 

Finally the poloidal variation of the potential is determined by quasineutrality, 

Z [ dh = '-^n^, (5) 



e 



giving e0?7Te ~ {B/Bp)6i. 

2.1.2. Turbulent distribution function and potential. The turbulent piece of the ion 
distribution function is obtained using the gyrokinetic equation (see Appendix |A]) 

Dfl 
Dt 



flvli - ■ 'VrH^i 



' [ui dip \T, 2) Ti dip BT, dip 
- ^ {uh + VM + vc) ■ VK(r ). + 5^ («b + Vm) ■ Vn{<P%, (6) 

where D/Dt = dt + -Rfi^C ' is the time derivative in the rotating frame, u = 



±^2[e — fiB + ~ ±a/2(£ — fiB) is the parallel velocity in the rotating frame, 

VAf = (/w/^j)b X Vr-B + {u^/Qi)h X (b • VrI)) are the VB and curvature drifts, 
Vc = i2un^/ni)h X [(Vi? X X b] is the Coriolis drift, v^^^ = -{c/B)Vn(pT x b and 
= -{c/B)'VR{(j)^^)i X b are the neoclassical and turbulent E x B drifts, C-plhf'} 
is the linearized ion-ion collision operator, and (...)« = {27i)~^ ^ dip {. . .)\-£i^E,fi,t is the 
gyroaverage holding the ion gyrokinetic variables R = r + fij~^v x h + . . E, fi and t 
fixed. The function that enters in the collision operator is 



' M \ Ti deo B dfio 

Here the subscript g in = /*'^(Rc,, w ^/2, v'j_/2B, t) indicates that we have replaced the 
variables R, e and /i by Rg = r + ^^~^v x b, t>^/2 and v']_/2B] similarly, the subscript 
in fMifl = /M^(^(r),t;V2,t) and H^l^ = H^^{iP{v),e{v),v^ /2,vl/2B,t) indicates that 
we have replaced the variables R, e and by r, and v']_/2B. 
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The equation for electrons is of the same form as the one for the ions, giving 



tb 



Dt 



+ ( Mb + vm + v^^ 



ptb 
' e 



= — V 



,tb 



1 dn^ 



Me 



H — 7^ [uh + Vm 



Tie dip \ 

) • VR(r) 



1^ 



Me 



(8) 



where vjvf = —{^,/VLe)h x VrB — [u^ /VL(,)h x (b • VrI)) are the V-B and curvature drifts 
for electrons, vg' = -(c/5)VR(0*^)e x b is the turbulent E x B drift, C'i] {hf} is the 
linearized electron-electron collision operator, cf^{h}'^^ hf'} is the linearized electron-ion 
coUision operator, and (. . .)e = (27r)~^ § dip{. . .)\^,E,ix,t is the gyroaverage holding the 
electron gyrokinetic variables R = r — ilj^v xh -\- . . E, ji and t fixed. The electron 
distribution function that enters in the collision operator is 

.^^tb _ ^^tb^^) 



leg "1 Tj, /Me,0- 



(9) 



The subscript g on ff^ = /*^(Rg, f ^/2, t'^/2i?, t) indicates that we have replaced the 
variables R, e and /i by R^, = r — f^^^v x b, and v'j_/2B; similarly, the subscript o 
on fiMefi = fMe{}pij)-,v'^/2, t) iudicatcs that we have replaced the variables R, e and ^ by 
r, f ^/2 and v']_/2B. If we were to neglect the effect of the trapped electrons, the solution 
to this equation for time and length scales typical of ion turbulence would simply be 
the adiabatic response ~ (e(0*'^)/Te)/Me because of the high parallel speed of the 
electrons. 

Finally, the electrostatic potential 0*^ is obtained from the quasineutrality equation 



M 



Ml 



Mi,0 



deo 
d\ft^g. 



B djiQ 



+ / d% 



fMe,0 

(10) 



2.1.3. Second order, long wavelength distribution function. The long wavelength pieces 
and il*2 are given by 



f2Me 



d^vS" 



+ 



- 1 



d'^vS' 



V 

where a — nc, tb. In the preceding. 



Me 



f. 



Mi 



rii 



(11) 



+ 



Me 



+ 



1 dpi 

Pidijj ' \ T, 

MIufMi dQc 
BTi dip 
ZedH^ 
'M~dr 



5\ i^m 

2) Tidijj 



Vm • VrV^ - Vm • VRi7°i 



«b-VR0r + Cf){if-,//-}, 



fMi (vc - -|VR0r X b) • VrV' 
ZefMi 



«b-VR0r + VM-VR0r) 

(12) 
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Table 4. Contributions to transport of momentum. 



n Size [{B/Bp)5^iPiR\Vi)\] Dependences 

11% (fip/fi)A„,(5ri for A„rf> d^n^, ^ud, d^T^, d^n^, d^T,, d^T,, d^rie 

1 for Aud<iB/Bp)6., 

K 1 d^T,,d^n,,d^T,,dlT,,dln,,dlT, 



U'^\ Aud{qRi^ii/vu)6~^ for Aud>{B/B.p)5i S^fi^ A«d, d^Ti, d^T^, d^Ue, d^Ti, d\ 
{B/Bp){qRvu/vu) for Aud<{B/Bp)5i 

{B/Bp) {qRvii/vu) 9^T„ d^n,, d.^T„ d^T,, din. 



with C}" the full bilinear ion-ion collision operator, and 

. _M V. . ( « <;.V^^>.) . f (.6 . .„) . v.(,.^)) J . ,13, 

2.2. Calculation of the momentum transport 

The radial transport of toroidal angular momentum {{RC' Pj ■V'?/')v)t is given in 
equation (39) of [18] that we reproduce here for convenience 



2Ze 
M^c^ 1 9 



1 rfV (a.{/.})T^'(v'-C)'^ 

This expression is derived in [18] . In Appendix [HI we present an alternative proof that 
makes clear the convenience of using form (fT4l) . 

Using that for B/Bp ^ 1, i?v ■ ^ ~ Iv\\/B, and employing the decomposition of 
the ion distribution function in subsection 12.11 we find that ([TH) can be rewritten as 



with 



m- P^ ■vv^),)t = n% + K + m + nr + (15) 



Lth 



dC J -"^ \ B 



Iv\\ 



^-i = Mcir-^ d^vf:^[-^ + Rn, ) , (16) 



^1 



T 



' T / ?/) 
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M^c I r . . 



d'vC^^H]^^^—^) , (17) 



2Ze 




^-1 = ( / d'v C^'{H^^ + Hll,}^ ) (18) 

and 

MV 1 d .J r ., „„,...„ 



(19) 



Recall that the subscript g indicates that R, e and fi have been replaced by = 
r + fi~^v X b, and vj_/2B, and the subscript o that they have been replaced by r, 
and v\/2B. In Table H] we summarize the size of all these contributions compared 
to the reference size {B / Bp)5^piR\Wip\, and we write what they depend on. To obtain 
these dependences, we use equations ([3]), (jl]), ([5]), ([6]), ([8]), ffTOj) and ffTTl) . The size 
estimates are taken from [18], where the turbulence was assumed to not scale strongly 
with Bp/ B. We use A^^ to denote a measure of the flux surface up-down asymmetry. It 
ranges from zero for perfect up-down symmetry to one for extreme asymmetry. Notice 
that for extreme up-down asymmetry, II^j^ and IV^i clearly dominate. 



3. Discussion 



We finish by showing how this new formalism gives a plausible model for intrinsic 
rotation. Until now, models have only considered the contribution 11^^, with ff^ 
and (f)^^ obtained by employing equations (EI) and (ITOil without the terms that contain 
H"^^ and 05*^. This is acceptable for RVL(^ ^ {B/Bp)6iVti or high up-down asymmetry 
A„rf > {B/Bp)6i. In these limits, ^^^(S^fif , fi^) ~ -z/*^9^^]^ - r*^^^ + H^f^. To 
obtain this last expression we have linearized around d^Q(^ = and = for 
RVLq/vu -C 1. Here z/**^ is the turbulent diffusivity, F**^ is the turbulent pinch of 
momentum and Ilj^'^ ~ ^udRPiR\^M is the value of 11*:^'^ at fi^ = and d^VLc, = 0, 
and is zero for perfect up-down asymmetry when equations ([6]), ([8]) and ffTOj) are 
solved without the terms that contain Hf^ and 05^"^ [24j. Notice then that imposing 

{{RC- Pi ■ViIj)^)t ^ n*^ = -z/*^9^fic - r^'^^c + njfd = O gives intrinsic rotation only 
for up-down asymmetry or if momentum is pinched into the core from the edge. 

The complete model described in this article includes contributions that have 
not been considered before. On the one hand, the gyrokinetic equations ([6]) and 
( ITOj) have new terms depending on Hf^ and 05^^, giving IV'\ ~ —i'^^d^Q(; — F**^!]^ + 
nj;^ -|- n^j^ Q, where FI^^ q ~ {B / Bp)5fpiR\V4'\ is a new contribution due to the new 
terms in the gyrokinetic equation. In addition, there are the new terms n*o^, Ul\ 
and IIq'^. As we did for II^j^, we can linearize U'^i^d^fl;^) around d^Q(^ = to find 
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-I- TT'^'^ 

{B/B,)\qRp,JvuW,V^R\V^|J\ 



AUB/Bp){qRuu/vu)SfpiR\Vi:\ and 
Combining all these results and imposing that 



where 11"^ 



Pi 



VV^)^)t 

■0(1 

dip' 



0, we obtain 



n 



int 



exp 




^tb _^ ^nc 



+n(lv.=V'a exp ^ j dip' 



■^tb 



^tb _^ ^nc 



(20) 



where V'a is the poloidal flux at the edge, r2^|^=^^ is the rotation velocity in the edge and 
nint ^ ]-[t,b ^ ntb^^^ ^ ntb ^ ^ nnc^^^ ]-[nc_ Notice that this equation gives a rotation 

profile that depends on 11'°* that in turn depends on the gradients of temperature 
and density, and the magnetic field geometry. The typical size of the rotation is 
~ iB/Bp)6iVti/R for Aud<iB/Bp)6i and fi^ ~ Am/R for A„d> (5/5^)5^. 
This new model for intrinsic rotation has been constructed such that the pinch 
and the up-down symmetry drive, discovered in the high flow ordering, are naturally 
included. By transforming to the frame rotating with we have made this property 
explicit. 
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Appendix A. Equation for the distribution function in the rotating frame 

In this Appendix we derive equations ([3]), (E]) and ( fTTl) for the different pieces of the 
ion distribution function, equations (jl]) and ([8]) for the different pieces of the electron 
distribution function, and equations ([5]) and ( ITOi) for the different pieces of the potential. 
These equations are valid in the frame rotating with angular velocity and we deduce 
them from the results in [18], obtained in the laboratory frame. 

In reference [H] we showed that in the limit Bp/B ^ 1, assuming that the 
turbulence does not scale strongly with Bp/B, the ion distribution function is given by 

/,(R,E,/i,t) = /A«(V^(R),i?,t) + i^r(^(R),^(R),^,/i,^) + i^2'=(V'(R),^(R),i^,/i,t) + 

Fl^iipiR), e{R), E, /i, t) + /f (R, E, /i, t), where the size of these different pieces is i^-^" ~ 
{B/Bp)6ifM^, i^T ~ {B/Bpf6UAH, ~ {B /Bp){vu/qRiyu)6UAH and /^"^ = /^^ + 
with fl}' ~ 6ifMi and ~ {B / Bp)6f fMi- The equations for the different pieces were 
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obtained from the gyrokinetic equation 



^ + R • Vr/. + = {Cu{m^, (A.l) 



where the time derivative R is 



R = M'b(R) + V,, - -|Vr(0). X b (A.2) 



and the time derivative E is 



E = -§[«'b(R) + v^,] ■ Vk(0).. (A.3) 



Here, u' = ±^y2{E — ^B) is the gyrokinetic parallel velocity in the laboratory frame, 
and 

= X VrS + ^b X (b ■ VRb) (A.4) 

are the V-B and curvature drifts in the laboratory frame. Equations (19) and (20) of 
[TS] for F^^ and equation (24) of [18J for F^^ are obtained from the long wavelength 
axisymmetric contributions to (lA.ip of order SifMiVu/d and {B/ Bp)Sf fMiVu/a, 
respectively. Equation (25) of \l8l for FI2 is also a long wavelength axisymmetric 
component of (lA.ip . In particular, it is the contribution of order SffMiVu/d that when 
the equation is orbit averaged does not vanish as uu — )■ 0. Equation (55) of [18] for 
fl^ is the sum of the short wavelength components of (lA.ip of order SifMiVu/ci and 
{B/Bp)8ffMiVti/a. 

In this article, we write the formulation in [18] in the frame rotating with velocity 
that is, we need to use the new gyrokinetic variable e = E — IQ(^u'/B. Thus, the 
new gyrokinetic equation is 

^ + R■Vnf^ + i^ = (A.5) 

The time derivative of the new gyrokinetic variable e is 

• de 

e = R-Vne + E—. (A.6) 



c f B"^ 

R = ub + — ^b + va, + vc - -Vr(0), X b + O TTT^^^ti ) , (A.7) 



In R, using u' = u + IQ^B, with u = ±^2{e - fiB + leads to 

B " ' ' 5 ' - 

with 

VM = X Vr5 + — b X (b ■ VRb) (A.8) 

the V-B and curvature drifts in the rotating frame, and vc = {^uIVLq/ BVti)\3 x (b- VrIo) 
the Coriolis drift. To obtain this expression for R we have used {u'Y = + 2IVLc_u/ B + 
0[{B/Bpf5fvl] to write v^^ = vm + + Ol^B"^ / Bl)5lvti]. The usual result for the 
Coriolis drift = (2'uf2^/i7j)b x [(Vr-R x (^) x b] can be recovered by realizing 
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that for Bp/B < 1, b = C + 0{Bp/B), b ■ VrI) = -V^R/R + 0[{Bp/B)R-^] and 
I/B = R + 0[{BVB^)R], giving 



^-b X [(VRi? X C) X b] - ^ 



bx (b-VRb) + 0(5^i, 



(A.9) 



f2j -Bfij 

In addition, using /b/E = RC+hxVip/B, <p = 0o+0r+02'+0*^> (0o)i = 0o(V'(R),^) + 
0(52Te/e), = (jyTW'R),e{R),t) + 0[{B/Bp)6fTje] and (^f-^), = O [(57^2)52^^^]^ 

we can simplify equation flA.7p to 

R = ^ib + RQ^C + VA/ + vc - |VR0r X b - -|Vr(0*^), X b + O {6S^) ■ (A.IO) 
The time derivative e in (1A.6|) can be written as 



e = E 



lu' on 



B 

To simphfy this equation we use 



lu' 

R ■ Vr^ - ■ Vr ( — 



Bu' 



E. 



{All) 



R - Vr 



lu' 

~b' 



u'h ■ Vr 



lu' 

~b' 



+ v'm- Vr 



lu' 

~b' 



B 



(VR(0)i X b) ■ Vr 



lu^ 

~b' 



With this result, obtained by using those that follow in ( 1A.14P and ( lA.lSp . and 
employing = 0^ + 0- + 0- + ^tb^ ^^^^^ ^ 0o(^(R),t) + OiSfTje), (0-), = 
0r(^(R), e{R),t) + 0[{B/Bp)6fTje], (05=), = 05^(^(R), ^(R), t) + 0[(EVi?2)54T,/e], 
= M + 0[{B / Bp)6iVti], we find 

' b(R) + VM + vc] • (VR0r + VR0r + Vr(0*").) 

/ ("2 3 



(vA/ - ^Vr(0*^), X b ) ■ Vr^ + O 



i ti 



(A.13) 



To obtain the result in ( ]A.12p . we have employed ■ Vr'?/' = u'h ■ V-r,{Iu' /^li); 



lu' 

^(VR(0).Xb)-VR( — 



B 



B 



Zel 
MBu' 
Zel 



b X VrS 



Zel 



/ib X VRi^ 



l\2 



-h X VrIu 



VR(0)i 



r\2 



-b X (b ■ VRb) 



Vr(0), + O( -^5.^;^^, 



Br, 



MBu 



(A.14) 



where we have used I/B = R + 0[{Bl/B^)R] and b- VrI) = -VrIu R + 0[{Bp/ B)R-^] 
and 



■ Vr 



Vr X (u'h) - u'hh ■ Vr X b 



R 



lu' 



B 



52 



'^5^vl (A.15) 



where we have used b ■ Vr x b ~ {Bp/B)a \ b ■ V^^Iu'/B) ~ {a/ R){Rvti/ qR) ~ 
{Bp/B)vu and Vr x (w'b) ■ Vn{Iu'/B) = Vr ■ [u'h x Vr(/mV5)] = Vr ■ [(Jm'/S) VrC x 
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Vn{Iu'/B)] + Vr ■ [{u'/B) ( VC X V^) X V^{Iu'/B)] = Vr ■ { VrC x V^[P{u'f /2B-'] } - 
d(^[{u'/R^B)Vip ■ V^{Iu'/B)] = 0. 

With equations ( 1A.5I) . ( lA.lOj) and ( 1A.13I) . we can now easily obtain equations ([3]), 
([6D and ^ for H^^, /f , and H^^. To obtain ([3]), we take the long wavelength 
axisymmetric contribution to (1A.5I) to order SifMiVu/a, giving 

uh ■ VRi/- + VM ■ Vr/m. - §^«b ■ VR0r = CiP{H^}. (A.16) 

This equation differs from equations (19) and (20) of [18], and gives a function Hf^ 
different from the function defined in |T8j. The reason is that fMi{4'(R),£) + Hf^ + 
Hf^ must be equal to the function fMi{ip{^),E) + FJ"^ + F^^ defined in [18] to the 
order of interest, but how the terms of first and second order in 5i are assigned to one 
or the other piece differs depending on the frame. For this reason, we have changed 
the name of the functions. The final result in ([3]) is obtained from (IA.16P by using 
VM ■ VrV' = Mb ■ V^{Iu/Vti) for u = ±^2{e - fiB + R^ny2) ~ ±y/2{e - fiB). 

Equation is the sum of the short wavelength contributions to (lA.Sp of order 
SifMiVti/a and {B/Bp)6f fMiVu/a. The equation is straightforward if we apply the same 
methodology as in |18j . 

Equation (ITT]) is found from the long wavelength axisymmetric components of (lA.SP 
to order SffMiVu/a,- Note that to this order we have the time derivative dtfMi [18]. Using 
^tfM^ = [nT^dtn, + {Me/Ti-?,/2)Tr^dtT,]fMi and realizing that dtn, = {Jd^vS''')^ and 
{3/2)dt{niTi) = (/ d^vS'^Me)^ + (/ d^vS'^'^Me)^, we find the final form in ([11]), 
and ( lT3l) . Here the integral (/ d^v S^^Me)^ gives both the divergence of the turbulent 
radial energy transport and the turbulent heating. Similarly, (J d^v S^'^Me)^ gives the 
divergence of the neoclassical radial flux of energy. The equations for Hf2 and are 
obtained in the same way as equations (24) and (25) in [18j, i.e., the equation for is 
the axisymmetric long wavelength component of (lA.SP of order {B / Bp)6f fMiVu/a, and 
the equation for Hj^ is the axisymmetric long wavelength component of order SffMiVu/o, 
that when it is orbit averaged does not vanish as ua — )■ 0. 

The equations (jlj) and (jH]) for the electron distribution function in the rotating frame 
are derived in the same way as the equations for the ion distribution function. The only 
differences are that the Coriolis drift vc and the term in ( ]A.13I) that is proportional to 
d^^li^ are small by ^JmjM and hence negligible, and that we include the electric field 
E"^ driven by the transformer, leading to a modified time derivative for the energy 

e = -±uh ■ E-^ + -[Mb(R) + vm] • (VR^r + VR(0**^)e) . (A.17) 
mm ' 

Finally, the equations for the different pieces of the electrostatic potential ([5]) and 

(ITU]) are easily deduced from the results in [IB] by realizing that moving to a rotating 

reference frame does not modify the quasineutrality equation. 
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Appendix B. Derivation of equation ( 1141) 



In this Appendix we derive equation witli a procedure different from tlie one 
employed in [T8j. Tliis new derivation sliows tlie connection witli calculations that 
split the off-diagonal components of the viscosity into gyroviscosity and perpendicular 
viscosity [25]. The derivation presented here and the one in [18] lead to identical results 
(as they should), but we believe that this new approach emphasizes the advantages of 
the formula in ( JT^ . 

We begin by using RC, = Ih/B — B^^h x Vip to write 



{{RC- P^ ■V^)^)t 



M 

'b 



dh' fiV^ ■ Vip Iv'u - (v' X b) • 



MI 



^ / d\' UViIj ■ vlvl - (v' X b)(v' X b) • (b X V^) 



(B.i; 



Since the transport of toroidal angular momentum needs to be known to order 
{B / Bp)5fpiR\Vip\, evaluating it directly from this equation would require knowing /j to 
order {B /Bp)6f fMi, and it is easy to see from our decomposition of the ion distribution 
function given in Section 12.11 that we cannot calculate the ion distribution function to 
that order. To circumvent this problem, we use exact moments of the Fokker-Planck 
equation to write the two integrals that appear in the transport of momentum as 

_M d 



M I d\' f,vW 



d'v'fXiVxh) 



V ■ ( / d\' /iv'v V 



Mc 
M 

a 



b ■ V(/) / dh' fiV X b + ( V(/) X b) / d^v' fiv'u 



V— I dVC4/,}t;[|(v'xb), (B.2) 
where we have used the x b) moment of the Fokker-Planck equation, and 

y j c^V/, [vlvl-(v'xb)(v'xb)" 

T\/T ^ V r ^ 

rfV/, (^vl(v'xb) + (v'xbK 

V ■ f / dh' /,v W 



M d 
M ^ 



I -bb 



^fl-bb). 

4a 

Mc 



X b 



4B 

Mc 



d\'h 



(V(/)xb)v^ + v'^(V0xb) 



45 



d^v' fi V±0(v' X b) + (v' X b)V_ 
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(B.3) 



where we have used the v^(v' x b) + (v' x b)v^ moment. These are the 
standard expressions used to calculate the perpendicular-parallel and perpendicular- 
perpendicular components of gyroviscosity and perpendicular viscosity [25]. 

Upon coarse grain averaging, the terms with time derivatives in (IB. 21) and ( IB. 31) 
are smaller than because dt becomes of order dfvu/a and the lowest order 
distribution function is a Maxwellian, making the velocity integrals over f ' | (v' x b) and 

v^(v' X b) + (v' X b)v^ vanish. As a result, the time derivative terms are negligible, 
giving 



MI 
/ McI 



(b X W) ■ 



d'v' (/.)tv'v'v' 



+ 



^-^b-V0 / rfV/,(v'xb)-VV^ 
/McI 

MI 



-(V0xb)-V^ / d'v'f.v'n 



and 



M 
2B 



v^v^ - K X b)(v' X b) ■ (la X W^jj) 



M 



+ 



\ABa 

I M 

Mc 



V- dh' (/.)tv'vv 



V'/ T 



h X Vil)) ■ 



V - ( / rfV (/.)tvVV 



252 



/ Mc 



(b X V^) 



\2B^j ^ ^'^^t^"^' ^) ■ ^^H(V0 X b) • V^] 
/ M 



>.). 



\ ABVLi 



v' ■ vi)Y - (v' X b) ■ v^A 



(B.4) 



(B.5) 



We could use these two equations to evaluate ((-RC' Pj ■V^/')^)t instead of the form in 
( |T4l) . If we did so, we would still need to evaluate the first term on the right side of ( 1B.4P 
and the first and second terms on the right side of (IB. 50 to order ^fA^dPi-Rl VV'I for an 
up-down asymmetric tokamak with A^^ ^ {B/Bp)6i, and to order {B / Bp)6fpiR\Vip\ 
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for an up-down symmetric tokamak with Aud<{B / Bp)Si. In the up-down asymmetric 
case, the dominant contribution to the first term on the right side of flB.4p and the first 
and second terms on the right side of ( 1B.5I) is due to Hf^, 

MI 



BVt, 



(b X V^) ■ 



MI 

I MI 

M^c 
2Z7 



) ■ Vb • (b X VV^) J d'v Hl^v\\{vi - v^{] 
(bx V^)-V(^y" d\ Hl^v\\vl^^ ^ 



VlnS 



(B.6) 



This term is of order {B/Bp)Aud5fpiR\V%lj\, suggesting that the transport due to the 
neoclassical piece Hf^ does not scale with coUisionality and that it is larger than the 
turbulent piece 11^]^ by a factor oi B / Bp 1 (see Table H]). In fact, this contribution 
cancels to this order with other terms in ( 1B.4P and (IB.Sp . as we will show in equations 
f lRTl) . flRSl) . f ira . f lRTOl) and flRTTll below. The final result of the cancellation is 
that the neoclassical pieces of the distribution function give a contribution of order 
{B/Bp){qRv,Jvu)Kd5tPiR\Vxl^\ for A„, > {B / Bp)5,. 

In the up-down symmetric case, similar problems appear. Obtaining the transport 
of toroidal angular momentum to order {B / Bp)5^PiR\Wip\ requires calculating the long 
wavelength piece of the ion distribution function {fi)^ to order {B/ Bp)6f Jmi for the 
first term on the right side of flB.4p and to order {B / Bp^df fui for the first and 
second terms on the right side of ( IB. 51) . In section 12.11 we show how to calculate 
Hf^ ~ {B/ BpYSffui and H}^ ~ {B / Bp){vu/ qRuii)Sf Jmi, so in principle, it is possible 
to evaluate these terms to the appropriate order. Note, however, that equations (]B.4p 
and ( IB. 51) suggest misleading scalings for the transport of momentum. For example, 
using 11^2 ~ {B / BpY5f fui in the first and second terms on the right side of (IB.5p . 
we obtain that a purely neoclassical piece can give momentum transport of order 
[B / Bp)5^PiR\Vilj\, that is, transport that does not scale with coUisionality. Similarly, 
integrating over if*2 ~ {B / Bp){vti/qRi'ii)Sf fui in the first term on the right side of (]B.4p 
gives momentum transport of order {B / Bp){vti/qRi'ii)6fpiR\V'ip\, i.e. it scales inversely 
with coUisionality. In fact, these contributions to the transport of momentum vanish by 
themselves or when combined with other terms in ( IB. 41) and ( 1B.5P as we will show in 
( ]B.7|) , ( IB.8|) , (]B.9|) , (]B.10|) and (]B.11|) . In equation (IT^ these cancellations have already 
been taken into account. There are advantages to this. For example, if we decide to 
use (IB. 41) and (IB. 5 1) instead of ( ITil) . 11^2 ~ {B / Bp){vti/ qRi'ii)5f fui must be calculated 
to order {B / Bp)5f fui < {B / Bp){vti/ qRi'ii)5f fui, that is, to more precision than is 



necessary. Equation ( IT^ . on the other hand, makes explicit that the first term in the 
right side of ( IB.4P vanishes when the lowest order piece of H^^ ~ {B / Bp){vti/ qRi'ii)5f fui 
is integrated over and combined with other terms in ( IB. 41) and (IB. 51) . 

In what follows we show how to obtain ( IT^ from (IB.4P and (IB.5p . Using V^/'V^/' = 



Sources of intrinsic rotation in the low flow ordering 



18 



|VV'P(I -bb)-(bxVV')(bxVV'),bxVV' = /b-i?SCand V(i?C) = (Vi?)C-C(Vi?), 
we obtain 
/ MI 



\bvl, 



■(b X W) 



V ■ / d'v' (/.)tvVV' 



■b 



/ M 



M^c 1 (9 



(b X W)) 



2Ze 1/' (9^ 
/ MP ^ 
\2BQi 

\ 4Bf]i 



;--\/'(i?MdV(/,)T(v'-C)V-VV' 



V • / d^^;' (/^)tvVv' 



I -bb 



ciV(/OTvVv 



/Mc/, 



h V(j) dh'fiiv' X b) • VV^ 



+ 



Mc 



t^V/4(v' X b) • Vip][{V(f) X b) • V^] 



Mc (i?2^-V0 / dV/i(v'-C) 



' /AfclW 
\ 2^2 



T 



and 

"MJ 



(iV {Cu{fi})^vl^{v'xh)-vA 



+ 



/ M 

M^c 
~2Z^ 



(v' • VV')^ - ((v' X b) • VV^ 



(B.7) 



(B.8) 
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+ 



\ 25^],• 



/ \2 



/ \2 



(B.9) 



Furthermore, we will use the {MP /2BVt.-^{v'^i^'^ and {M\Vi)\^ / ABVLi){v'^_Y moments of 
the Fokker-Planck equation to remove most of the J d?v' (/i)Tv'v'v' moments, 



/ MP 



\2Ba 



b- 



V- (^y f/v(/,)TvVv' 

/ MP 



■b 



/McP, 



fi2 



b- V0 / d'v'fA 



+ 



\2Bn 

d I MP 



dt \ 2BVLi 



/ 



\ 25a 



9pi 

9t 



(B.IO) 



and 



45a 



+ 



(B.ll) 



Combining f lRT]) . f lRSj) and flR9|) . and then using fIB.lOp and fIB.lip . we finally obtain 



{{RC- Pi ■VV^)^)t 



Mc 
2Z'e 



dpi M^c 1 d 



dt 2ZeV'dip 



V'(R' I rfV(/.)T(v'-C)V-V^ 



+Mc (R\-V<P / dV/,(v'-0 



M^c 
2Z7 



(B.12) 



To obtain the final result in (III]) , we just need to rewrite the integral {R^ J d^v' (/i)T(v'- 
C)^v' ■ Vip)ip using the Fokker-Planck equation as is done in [18] . 

Finally, note that in writing (HM we have used the expressions ( ]B.2p and ( IB.SP 
that give the gyroviscosity and perpendicular viscosity, and ( IB.lOp and ( 1B.11|) that are 
equations for the parallel and perpendicular pressure. It is necessary to combine all 
these equations to explicitly show the cancellations mentioned below equation (IB.Sp . 
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